function IR=impkw(M,PI,G,ylist,xlist,lpd,shock,nir);

% Matlab File: impkw.m
% calculation of impulse response (for nir periods) using the Markov
% decision rules. The inputs and outputs are summarized by the statement:
%
%    IR=impkw(M,PI,G,ylist,xlist,lpd,shock,nir) 
%
% where G is the matrix in the driving process specification, 
% x(t) = Q delta(t)  and  delta(t) = RHO * delta(t-1) + G e(t), 
% lpd indicates the location of the predetermined variables and 
% shock indicates the value that e(t) should take. If no shock vector
% is specified, then a menu allows a choice of shocks.

nM=rows(M);
nPI=rows(PI);

if (nargin<7)
   clc
   disp('No shock vector specified. User can supply shock if G=eye(M)')
   disp(['The matrix M contains ', num2str(nM), ' elements'])
   klist=ylist(lpd,:)
   nk=rows(klist);
   for i=1:1:nk;
      sklist(i,:)=[sprintf('%2.0f',i), ' ', ylist(lpd(i),:)];
   end
   disp(['The first ', num2str(nk),' elements are predetermined variables'])
   disp(sklist)
   nx=rows(xlist);
   for i=1:1:nx;
      sxlist(i,:)=[sprintf('%2.0f',i+nk), ' ', xlist(i,:)];
   end
   disp(['The next ', num2str(nx),' elements are exogenous variables (x)'])
   disp(sxlist)
   shk=input('select shock # ');
   szshk=input('choose size of shock ');
   shock=zeros(1,nM);
   shock(:,shk)=szshk
   pause
   nir=input('choose length of impulse response ');
end

   
nM=rows(M);
nPI=rows(PI);

% IMPULSE RESPONSES

% The calculation of impulse reponses involves two elements.  First, it is 
% necessary to specify initial conditions for the controlled state(k). Second,
% one needs to specify an values of the exogenous variables, i.e.,
% an (nM x 1) vector of initial conditions.

if (rows(shock)<cols(shock))
 disp('transposing shock vector')
 shock=shock';
end

if max(size(shock))~=max(size(M));
 disp('inadmissable run')
 IR=0;	
else

S = shock;

% Define initial value of state vector for use in later work,
% since S will be updated in computation of impulse response.

SI=S;



% Initialize the impulse response matrix

IR=zeros(nir,nPI+1);

% Calculate the impulse reponse for the nir periods.

for i=1:nir;
  %Index of Impulse Response Period;
    IR(i,1)=i;
  %Observables;
    IR(i,2:nPI+1)=(PI*S)';
  % Update the state vector:;
    S = M*S;
end

end